PolyOptimalIRT:

library(mirt)
library(psych)

titlestr  <- "hads-7-anxiety-screen"

U <- scan("data/hads.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

#  drop cases containing missing data
keepindex <- U[,14] < 4
N = sum(keepindex)
#  N = 807

U <- U[keepindex,]
# gives names for mirt
colnames(U) <- paste("I", seq(1, ncol(U), 1), sep = "")

scree(U)

hadsGrm = mirt(data = U, model=1, itemtype = 'graded', SE=T)
hadsGrm2 = mirt(data = U, model=2, itemtype = 'graded', SE=T)
summary(hadsGrm2, suppress=0.25)

anova(hadsGrm, hadsGrm2) # p = 0 indicates better fit for 2 dimS
# save(U, hadsGrm, file="data/hads_fittedmodel_mirt.RData")
for(i in 1:ncol(U)){
    ItemPlot <- itemfit(hadsGrm,
                        group.bins=15,
                        empirical.plot = i,
                        empirical.CI = .95,
                        method = 'ML')
    print(ItemPlot)
}
itemfit(hadsGrm)

# confirmatory mirt -------------------------------------------------------
fstruct <- "F1 = 1,3,5,7,9,11,13
F2 = 2,4,6,8,10,12,14"
hadsGrmConf = mirt(data = U, model=fstruct, itemtype = 'graded', SE=T)
summary(hadsGrmConf)




# Random ------------------------------------------------------------------

grmpars <- coef(hadsGrm, simplify=T)$items
grmpars <- grmpars[, c(1, grep("d[1-9]", colnames(grmpars)))]
sampleGRM <- simdata(grmpars[, 1], grmpars[, 2:ncol(grmpars)],
                     itemtype = 'graded', Theta=matrix(c(0, -2, 2), ncol = 1))
fscores(hadsGrm, method = "EAP", response.pattern = sampleGRM)[, "F1"]

plot(hadsGrm, type = 'trace', which.items = 1:4, facet_items=T, theta_lim=c(-6, 6))
plot(hadsGrm, type = 'trace', which.items = 1:4, facet_items=T, theta_lim=c(-3, 3))
plot(hadsGrm, type = 'infotrace', which.items = 1:4, facet_items=T, theta_lim=c(-4, 4))
plot(hadsGrm, type = 'score', facet_items=T, theta_lim=c(-4, 4))